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We study in experiment and with computer simulation the free energy and the kinetics of vacancy 
and interstitial defects in two-dimensional dipolar crystals. The defects appear in different local 
topologies which we characterize by their point group symmetry; C n is the n-fold cyclic group 
and D n is the dihedral group, including reflections. The frequency of different local topologies 
is not determined by their almost degenerate energies but dominated by entropy for symmetric 
configurations. The kinetics of the defects is fully reproduced by a master equation in a multi-state 
Markov model. In this model, the system is described by the state of the defect and the time 
evolution is given by transitions occurring with particular rates. These transition rate constants are 
extracted from experiments and simulations using an optimisation procedure. The good agreement 
between experiment, simulation and master equation thus provides evidence for the accuracy of the 
model. 


I. INTRODUCTION 


The microscopic dynamics and interaction of defects, 
like dislocations, vacancies and interstitials are key to 
a variety of macroscopic phenomena of materials mia. 
In two dimensional systems [SHU the peculiar two step 
melting is a result of dislocation and disclination inter¬ 
actions and melting is mediated by the formation and 
subsequent dissociation of dislocation pairs [MS]. Dis¬ 
locations pairs in 2d-crystals also form as the result of 
spontaneous clustering of interstitals and vacancies intro¬ 
duced into the systems PI- The dynamics and interac¬ 
tion of interstitials and vacancies is even related to exotic 
phases such as supersolidity mm- Nevertheless, even 
though individual defects and their displacement fields 
are well described by elasticity theory, the precise kinet¬ 
ics and the sign of the interaction cannot be captured by 
the theory ns 120] . The open question is how the de¬ 
fect kinetics emerges from the non-linear effects near the 
defect centers. Video microscopy of two dimensional col¬ 
loidal crystals together with optical tweezers now allow 
to investigate such nonlinear effects with single particle 
resolution in real time mm- 

Here, we show that the dynamics of defects in two- 
dimensional dipolar crystals can be fully described by 
a sequence of jumps between states which are defined 
by the local displacements in the vicinity of the defect 
centers. The equilibrium probabilities of the states (i.e. 
the populations) are a result of the interplay between 
entropic and energetic contributions. We find that the 
different contributions can be understood quantitatively 
from statistical mechanics by a harmonic expansion of 
the energy around the minima corresponding to the var¬ 


ious defect states. The kinetics of the defects follows a 
master equation for which we measure the transition ma¬ 
trix from a long experimental trajectory with the aid of 
an optimization routine. The results from experiments 
are compared with the results from Monte Carlo simula¬ 
tions. The equilibrium probabilities from simulation and 
experiment are in qualitative good agreement but show 
systematic differences for stiff crystals. 
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FIG. 1: Typical snapshots of the colloidal crystal obtained 
from the experiment containing an interstitial (upper row) 
and a vacancy (lower row), respectively. The motion of the 
defects consists of a sequence of transitions between states 
with different symmetries given by the dihedral D n and cyclic 
C n point groups in 2D. The upper left index denotes intersti¬ 
tial versus vacancies and the upper right counts the numbers 
of seven-fold particles (in case of interstitials) and five-fold 
particles (in case of vacancies). Particles with 4, 5, 7, and 8 
neighbors are shown in blue, orange, green, and purple, re¬ 
spectively, and particles with 6 neighbors are shown in black. 
Note that the symmetry of the most frequent low tempera¬ 
ture configurations (second column) differ for interstitials and 
vacancies. The vacancies show only cyclic symmetry for low 
temperatures. 
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II. EXPERIMENTAL SETUP AND MODEL 


We consider a system of super-paramagnetic colloidal 
particles (4.5/rm in diameter) confined by gravity to two 
dimensions on a flat interface. Two situations were re¬ 
alized, a flat water/air interface in hanging droplet ge¬ 
ometry and a solid interface of a glass substrate. The 
only difference in both realizations is a slightly enhanced 
self diffusion coefficient of colloidal particles at the wa¬ 
ter/air interface on the expense of an increased equili¬ 
bration time after mounting the sample compared to the 
solid substrate. The latter is due to a nontrivial regu¬ 
lation of the curvature of the droplet [26;. Video mi¬ 
croscopy and digital image analysis provides the position 
of the particles at a frame rate of about 1 Hz which is 
fast compared to the Brownian timescale of 50 sec. The 
colloidal ensemble is described by the Hamiltonian 
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TABLE I: Classification of interstitial (I) and vacancy (V) 
states based on the number of particles with 4, 5, 7 and 8 
neighbors in the vicinity of the defect. 


III. DEFECT CLASSIFICATION AND 
PATTERNS 
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with the pair interaction for particles at distance r 


f3V(r) = 


3 3 / 4 r 

(2tt) 3 / 2 



( 2 ) 


Here, /3 = 1 /ksT with the Boltzmann constant ks and 
temperature T. The distance is given in units of a, the av¬ 
erage inter-particle distance of the triangular lattice. The 
dimension-less parameter T = /3(/io/47r)(xiL) 2 (7rp) 3 / 2 
defines the phase behavior of the system and can be in¬ 
terpreted as an inverse temperature, tuned with the mag¬ 
netic field H . Here, po is the permeability of vacuum, x 
is the magnetic susceptibility of the particles and p is the 
2d particle density. The experimental setup is described 
in detail in Refs. 16, [24, 126] , On the solid substrate a 
vacancy is prepared by trapping a colloid with an opti¬ 
cal tweezer (20mW, lOOx tweezer-objective from above 
the sample, NA 0.73, Ar + -Laser) and pulling it along¬ 
side a lattice line out of the field of view (at least more 
than 15 lattice spacings). Correspondingly, an intersti¬ 
tial is created by pulling a particle from the far field into 
the center of the field of view of an otherwise defect free 
crystal. Samples confined at the water/air interface of¬ 
fers additionally the possibility to shoot particles out of 
the interface by light pressure with a strong laser pulse 
(500mW). Since the sound velocity is of the order of a 
few mm/s the distorted lattice relaxes rapidly and mea¬ 
surements are started after about a minute. Positional 
data were taken at three different interactions strength 
T = 120,140,156 with more than 3,000 configurations in 
the case of T = 120 and 156 and more than 20,000 con¬ 
figurations in the case of T = 140 for both interstitials 
and vacancies. For T = 140 the creation of the defect 
was repeated frequently, annealing the crystal and equi¬ 
librating the systems in between for several days up to a 
few weeks. 


In equilibrium, the vacancies and interstitials exist in 
different, almost degenerate states with different topolo¬ 
gies EI1I1H]- The most frequently appearing vacancy and 
interstitial states are shown in Fig. [l]and we classify the 
defects according to their point groups in two dimension. 
C n denotes the cyclic group with n-fold rotational sym¬ 
metry whereas D n is the dihedral group with additionally 
n mirror axes. The upper left index denotes vacancies or 
interstitials and the upper right index counts the number 
of dislocations involved in the defect. In 2D a dislocation 
is a pair of a five- (colored orange in Fig. [l]) and seven-fold 
(colored green) coordinated particles characterized by a 
burgers-vector. The neighbor numbers in the vicinity of 
the defect centers are determined using a Voronoi con¬ 
struction |29]. The dissociation of dislocations is known 
to drive the melting transition in 2D [8]. The most fre¬ 
quent defect configurations are listed in table [T| and the 
index i is introduced to label the defects in the formula 
below and for the computations of transition rates. Con¬ 
figurations with larger numbers of dislocations exist and 
are numbered with others (i = 0). Note, that the num¬ 
ber of such defects vanishes for large T (see dashed lines 
in Fig. [6|. 

Interestingly, the typical symmetry of the configura¬ 
tion consisting of two dislocations (second column in 
Fig. 0 is not the same for interstitials and vacancies. To 
quantify the relative probability of dihedral and cyclic 
configurations we introduce the following procedure. We 
connect the five-fold particles with a line and the seven¬ 
fold ones too and measure the angle a between the two 
axis. For a = 7 t/2±7t/8 the pattern is classified as rhom¬ 
bic ( D 2 symmetry) and otherwise classified as cyclic ( C 2 
symmetry). Fig. [2] shows the probability density as heat 
map for vacancies (a) and interstitials (b). Red corre¬ 
sponds to large probabilities and green to low ones. For 
T > 120 the distribution is sharply peaked at about 
a « 7r/2 for the vacancies and the probability of the 
D 2 symmetry tends to zero. The ratio for D 2 symme¬ 
try increases to 30% at T = 80. Having this in mind 
we keep the notation v C\ for vacancies consisting of two 
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FIG. 2: Probability density of the angle a between the 5- 
5 and 7-7-bonds for vacancies (left) and interstitials (right) 
consisting of two dislocations for different values of T. Red 
indicates high probabilities. For a — 7 t /2 the symmetry is 
dihedral D 2 . The vacancies are peaked at a « 90° ± 38° for 
low temperatures (large T) with cyclic symmetry C 2 - 

dislocations in the following plots for clarity. For inter¬ 
stitials Fig.[2jb) the distribution is peaked at a ~ 7 t /2 
corresponding to dihedral symmetry D 2 for all interac¬ 
tion strength. Nonetheless the distributions widens for 
lower temperatures which might indicate a separate con¬ 
figuration C 2 . Using the same cutoff as for vacancies the 
ratio varies form 50% : 50% at Y = 80 to 75% 1 D 2 and 
25% t C 2 at T = 170. 

To visualize the configurations with respect to the un¬ 
derlying lattice of the crystal, Figures [3] and [4] show the 
probabilities of positions of particles with respect to the 
defect centers. All patterns reflect the symmetry of the 
crystal but the relative positions of particles with 4, 5, 
and 7 neighbors for the interstitials are arranged in par¬ 
ticular patterns shown in Fig. [3] and the 5, 6, and 8 
neighbored particles of the vacancies are shown in Fig. [4] 
Remarkably, the positions of the 5-coordinated particles 
in an interstitial defect and the positions of 7-coordinated 
particles are all degenerate within the 12 points in a star- 
like pattern. 

IV. EQUILIBRIUM DEFECT POPULATIONS 



FIG. 3: Relative positions of 4, 5, and 7-coordinated parti¬ 
cles in the vicinity of an interstitial defect with respect to 
the defect position on the underlying triangular lattice from 
computer simulations. The configurations (left column) 1 D 2 
(a), 1 D 2 (b) and 1 (c) consist of at least two 5-coordinated 

particles (green). The 1 D 2 state consists of an additional 
4-coordinated particle (yellow) while 1 D\ and 1 D\ state con¬ 
sist of additional 7-coordinated particles (orange). The rel¬ 
ative positions of 5-coordinated (middle column) and 4- and 
7-coordinated particles (right column) are measured with re¬ 
spect to the closest lattice point of the defect center and de¬ 
picted as a heat map. Red indicates large probability and 
white low probabilities. Clearly, the 5-coordinated parti¬ 
cles occupy always the same 12 lattice points with different 
weights depending on the state. The 5 and 7-coordinated 
particles also follow different patterns. 


The probabilities are normalized to JT P i = l for va¬ 
cancies and interstitials separately. The population Pi 
equals the fraction of time the defects spends in state i 
in equilibrium. This probability is estimated from a finite 
number of configurations, 


Following a trajectory x(t) of the system in experiment 
one can use the classification above to identify a trajec¬ 
tory of states. Here, x(t) denotes the configuration of the 
system, including the positions of all particles, at time t. 
In computer simulations, x(t) is the sequence of Monte 
Carlo configurations. For illustration, a typical trajec¬ 
tory taken from computer simulations is depicted in Fig. 
[5] From the trajectory, we determine the probability Pi 
of finding the defect in state i with 

Pi = {hi[x(t)}) = lim - f dt'hi[x(t')\. (3) 

T^oo T J 0 

The indicator function fo(x) is defined as 

1 if x is in state i, / 

(4) 

0 otherwise. 


1 x \ 

p m = h = — ^2 h m [x(iAt)]. (5) 

i= 1 

The configurations are sampled in experiment and sim¬ 
ulation at constant time intervals. These configurations 
may therefore be correlated and the error from this esti¬ 
mation, including correlated events is 

a 2 = ((h - (h)f) = T J2(Sh(i)6h(j)), (6) 

ij 

where (h) is the ensemble average of h, and 
5h(i) = h[x(iAt)\ — (h). The correlation function 
(■ Sh(i)Sh(j )) can be rewritten as 
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FIG. 4: Relative positions of 5, 7, and 8-coordinated parti¬ 
cles in the vicinity of a vacancy defect with respect to the 
defect position on the underlying triangular lattice from com¬ 
puter simulations as described in the caption of Fig. [3] A 
vacancy in state v D® (a), v C\ (b), or v D% (c) consists of at 
least two 7-coordinated particles (orange) and either a single 
8 -coordinated (violet) particle (a) or additional 5-coordinated 
particles (green) (b) and (c). In all states, the 7-coordinated 
particles are found at the same 12 positions (right) while the 
8 and 5-coordinated particles are found at various other po¬ 
sitions. 


H(\i-j\) = (6h(i)6h(j)) = (h(i)h(j)) — (h) 2 = (7) 

= C m ,m(Ati-j) ~ ( h ) 2 . 

Here, C m m are the auto-correlation functions as depicted 
in Figs. |8] - [Tl] Note, that in equilibrium H(i — j) = 



t 


FIG. 5: Typical trajectory of the state as a function of time (in 
units of Monte Carlo time steps) of an interstitial at V — 160. 
The state is identified in each time step from the neighbor 
numbers of particles in the vicinity of the defect according to 
Table I. 


H(\i — j |) and H{ 0) = (h}( 1 — (h)). Inserting Eq. ([7]) 
into Eq. © with t = \i — j | we find 


/TV - 1 \ N-l 

ff2 = y 2 E ff W- ff (°) -ysE^)- (8) 

\ t=0 J t =1 

In the limit of large TV, the last term, which scales 
quadratically and the constant term H( 0) can be ne¬ 
glected and the error can be estimated from 




A ff (°) 


»<i - m K mp_ (9) 


N 




Note that the term Ml = a 2 - 1S error f r0 m N 

uncorrelated measurements. In a general trajectory cor¬ 
relations increase the error by 2 77(5) * This factor 

can be associated with the correlation time t c in the sys¬ 
tem by 


N—l 


E 


m 

H( 0) 


1 1 

AtH( 0) 



tc_ 

A t 


( 10 ) 


Here, the correlation time is t c = / 0 °° H(t)dt and 

At is the sampling interval. Combining all this, the error 
from correlated trajectories can be written as the uncor¬ 
related error multiplied by the correlation time 


2 _ 2 
a — a o 


2 tc 
At' 


( 11 ) 


In the experiments, the trajectory is sampled every 
At = 0.92 seconds. The number of measurements dif¬ 
fers considerably for the various T for vacancy and in¬ 
terstitials. In particular, for interstitials 7Vp =120 = 670 
and TVp =140 = 31, 500 and for vacancies TVp^ 120 = 4, 950, 
^rCi40 = 23, 700 and TVp^ 156 = 9,900. The correlation 
time for different species ranges from t c = 20 to t c = 40. 

In the computer simulations, the trajectory is a se¬ 
quence taken from Monte Carlo updates. The number 
of particles is N p = 26 x 30 + 1 = 781 for interstitials 
and N p = 779 for vacancies. Interactions are cut off 
at a radius r c = 8ao, where ao is the lattice spacing. 
The displacement in each Monte Carlo step are chosen 
such that the average acceptance rate is approximately 
0.5. Sequences are sampled every 800, 000 Monte Carlo 
steps which corresponds approximately Is in real time. 
For all parameters, the total number of measurements is 
N = 750,000. 

Populations of interstitial and vacancy states obtained 
from computer simulations for T ranging from T = 100 
to 170 are shown in Fig. [6] The defect populations de¬ 
termined experimentally for interstitials and vacancies at 
T = 120,140 and 156 are in good agreement. For the va¬ 
cancies the small overpopulation for T < 140 of y Cf-type 
with the lowest symmetry (red squares) is attributed to a 
tiny shear within the sample since vacancies were mainly 
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created at the water/air interface. This interface is less 
stable compared to the solid substrate but vacancies can 
be created with less perturbations by pushing particles 
out of the interface using the light pressure of the laser 
pulse. 
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FIG. 6: Equilibrium populations of (a) interstitial and (b) 
vacancy states identified based on Table I as a function of 
r from experiments (filled symbols) and simulations (solid 
lines). The dotted line denotes the probability to find the 
defect in a state not given in the table which are negligi¬ 
ble for higher interaction strength T. Error bars depict the 
statistical error from finite numbers of measurements taking 
into account the correlations between the individual configu¬ 
rations. 


V. FREE ENERGY OF DEFECT STATES 

The probability of finding the defect in state i is 
related to the free energy of the state, which con¬ 
sists of an energy and an entropy contribution. For 
convenience we introduce the reduced energy W(x) = 

2 * 3/2 Ei<j ( a / r ij) 3 i which is related to the potential en¬ 
ergy V(x) by /3V(x) = TW(x). The probability density 
in configuration space for a given temperature and vol¬ 
ume is given by 


p(x) = Z~ 1 e~ rw ^ (12) 


with the partition function Z = f e~ r " (x) dx. The prob- 
ability of finding the defect in state i can then be written 
as ensemble average, 

= f e~ rw Whi(x) dx _ Zi —rAFi 
1 fe- rw ^)dx Z 

where Zi = f e~ rw( ^hi(x) dx is the partition function 
restricted to state i. Thus, A F{ = F{ — F is the difference 
between the free energy Fi = —T _1 In Zi of configuration 
i and the free energy F = —T _1 InZ of the system. Ac¬ 
cording to basic statistical mechanics, the free energy Fi 
is the sum of an energetic and an entropic contribution, 

Fi = (W)t - (14) 

Here, (W)i is the average energy given that the system 
is in state i, and Si is the entropy of state i. 

For sufficiently large values of T, at which configura¬ 
tions belonging to state i can be viewed as small fluctu¬ 
ations about a local energy minimum, the energetic and 
entropic contributions to the free energy can be com¬ 
puted analytically. In this regime, the energy of each 
state i is approximated as quadratic function centred at 
the minimum energy configuration x^, 

Wi(x) = TT(x-) + iujPiUi, (15) 

where = x — x° is the displacement from the minimum 
and T>i is the matrix of second derivatives of W evaluated 
at x?. The partition function of state i is then given 

by Zi = e~ rw ( x °) ^ ? where n is the number of 
degrees of freedom. From this expression it follows that 
the free energy difference A F k i = F k — Fi between two 
states k and /, which determines the relative population 
Pk/Pi, can be expressed as 


AF kl = 


~Iln§ = 

r Pi 


(16) 


xx7-/ o\ xxx/ 0\ det T) k 

= W(x k .) - W(xj) - -- In -——- 
v kJ v l) T 2 det Vi 

= AW k , - IA Ski. 

Thus, the free energy difference A F k i = F k — Fi de¬ 
pends linearly on 1/T. The intercept of this function 
with the y -axis then yields the energy difference AW k i = 
FF(x^) — Wi(x®) between states k and l and the slope 
equals the entropy difference A Ski- As shown in Fig. [7| 
this linear behavior is indeed observed in our simulations. 
Figure [7] also shows the energy and entropy differences 
obtained by linear fits of Eq. (16) to the simulation re¬ 
sults. As can be inferred from the figure the symmetric 
defect configurations v D® and 1 D\ have a positive en¬ 
tropy with respect to the respective lowest energy states, 
leading to a negative slope of the free energy vs. 1/T 
curves shown as black lines in Fig. [7] (a) and (b). Note, 
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that while the energy difference between states is small 
for all T, the entropy vs. energy ratio may change dra¬ 
matically as a function of T. This positive entropy dif¬ 
ference for the symmetric defect configurations causes an 
inversion of the population order for lower values of T. 
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FIG. 7: Free energy difference A Fki as a function of l/r for 
vacancies (a) and interstitials (b). The dashed lines are linear 
fits to the data obtained from computer simulations. Accord¬ 
ing to Eq. (16), the intercepts of these functions corresponds 
to energy differences A Wki shown in panel (c) for interstitials 
(black) and vacancies (red). The energies are given with re¬ 
spect to the lowest energy state which is v C% and 1 for va¬ 
cancies and interstitials, respectively. The slopes of the lines 
in panels (a) and (b) yield the entropy differences A Ski shown 
in panel (d) for interstitials (black) and vacancies (red). 


VI. DEFECT KINETICS 

On a coarse level, neglecting microscopic details, the 
motion of a defect can be viewed as a sequence of transi¬ 
tions between discrete states. This type of dynamics can 
be mapped onto a Markov state model governed by the 
master equation 

[K^t) - KjiPM. (i7) 

j/* 

Here, Pi(t) is the probability of finding the defect in state 
i at time t and is the rate constant for transitions 
from state j to state i. The general solution of the master 
equation can be written in matrix and vector notation as 

P(£) =exp[Kt]P(0). (18) 

Here, P(£) is the vector of probabilities Pi at time t and 
K is the matrix of transition rate constants Kij. While 
for short times P (t) depends on the initial conditions 
P(0), for long times P (t) converges to the vector P eq of 


equilibrium populations independent of time and initial 
conditions, 


P eq = lim exp[K£]P(0), (19) 

t—yoo 

For the matrix of rate constants, the condition of detailed 
balance holds with respect to the equilibrium distribution 

peq 

1 i 


PT = Kji 

PP Kij- 


( 20 ) 


In addition, the conservation of total probability requires 


Ku = — Kji. (21) 

J#* 


The transition rate constants can be calculated from 
a trajectory of states with the following procedure. To 
characterise the time evolution of the system, we intro¬ 
duce the correlation functions 


Cij{t) = P(i,t\j,0) 


{h i [x(t)]h j x.( 0)]) 
(/ij[x(0)]) 


( 22 ) 


The correlation function Cij(t ) is the conditional proba¬ 
bility of finding the defect in state i at time £, given that 
it was in state j at time t' = 0. The equilibrium proba¬ 
bility P® q of state i is the large time limit of P(i, t\j, 0), 

P^= lim P(i,t\j,0). (23) 


To obtain this equation, we have used the fact that 
for long times the state of the system at time t is sta¬ 
tistically independent of the state at time t — 0, i.e., 
(hi(t)hj( 0)) = (hi(t))(hj( 0)). The correlation functions 
Cij(t ) can be easily determined from trajectories ob¬ 
tained in experiments or simulations. This correlation 
function is then compared to the result of the master 
equation given the matrix of rate constants with 


Cij(t) = (exp[K£]P) mt ).. (24) 


Here, the initial vector P^ mt has a 1 in component j while 
all other components have a value of 0. This particular 
choice of initial condition implies that system is initially 
in state j with probability 1 as required by the definition 
of the conditions probability Cij. 

To determine the matrix of rate constants Kij govern¬ 
ing the dynamics of defects, we carry out an optimization 
procedure that minimises the difference between the time 
correlation functions measured in experiments or simula¬ 
tions and those predicted from the solution of the master 
equation. The target function of this optimisation is de¬ 
fined as 


B < K ) = EE[ C i:j (lAt) - C i:j (lAt, K) 

ij l 


( 25 ) 
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FIG. 8: Time correlation functions Cij(t) for interstitials ob¬ 
tained from experiment (solid lines) in comparison to the 
result of the optimization procedure of Kij (dashed) for 
r = 140. The symbols in the top right corner of each panel 
indicate the initial state j and colors indicate the final state 
i, where black corresponds to 1 D®, red to 1 green to 1 
and blue to v /|. For instance, the black lines in the top left 
panel represent the probability to find the defect in the 1 D® 
configuration at time t given that it was in 1 D® at time t — 0, 
while the red lines represent the probability to find it in 1 D\ 
at time t. 




FIG. 9: Correlation functions Cij(t) for interstitials obtained 
from simulation (solid lines) and from the optimization pro¬ 
cedure of Kij (dashed). Parameters and colors as in Fig. [8] 


Here, Cij(At) denotes the matrix of correlation functions 
determined from experiment or simulation. The argu¬ 
ment K in the target function E( K) and in the matrix 
of predicted correlation functions C^-(At,K) emphasise 
their dependence on the matrix of rate constants. The 
matrix of rate constants is optimized to best reproduce 
the observed time correlation functions. In this optimiza- 


FIG. 10: Correlation functions Cij(t) for vacancies obtained 
from experiment (solid line) and from the optimization pro¬ 
cedure of Kij (dashed). Parameters and colors as in Fig. [8] 



t 


t 


FIG. 11: Correlation functions Cij{t) for vacancies obtained 
from simulation (solid line) and from the optimization proce¬ 
dure of Kij (dashed). Parameters and colors as in Fig. [8] 


tion procedure the target function of Eq. (25) is then 
minimized iteratively. In each step of the iteration, a 
matrix element K^ is chosen at random and changed by 
adding a random amount A k. To satisfy the constraint 
Eq. (21 ) A k is also subtracted from matrix element Ka . 
The step is accepted if the target function has decreased. 
The iteration is stopped when the target function has not 
decreased for a certain number of steps. The optimisation 
procedure is initialized with the time derivatives of the 
time correlation functions Cij(t ) evaluated at t = 0, with 
Kjj ST = dCij(t)/dt | t =o- These initial transition rate con¬ 
stants correspond to the transition state theory estimates 
obtained for the dividing surfaces defined implicitly by 
the state classification introduced earlier. In the transi- 
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FIG. 12: Color coded optimized rate matrix K ij for an inter¬ 
stitial (a) and (c) and a vacancy (b) and (d) from a trajectory 
from experiment (top) and simulation (bottom) at T = 140. 


tion state theory approximation, correlated crossings of 
the dividing surface are neglected. Note, that the opti¬ 
mization procedure described above takes these transient 
short-time correlations correctly into account. 

Transition rate constants obtained with this opti¬ 
mization procedure for vacancies and interstitials are 
summarized in Fig. [12] The results from computer 
simulations and experiments are in very good agreement. 
The results show that the transition rates are not homo¬ 
geneous. The most dominant transitions are 1 D® — D\ 
for interstitials and v D\ ~^ v for vacancies. We also 
identify several transitions with a vanishing rate (e.g. 
v I f D 2 or v D 2 ~^ v Cf. We plan to investigate the 

exact microscopic mechanism for these rare transitions 
in the future. 


accurately determined from a second order expansion of 
the energy with respect to displacements. 

The kinetics of the defects is well described by a mas¬ 
ter equation in a multi-state Markov model. The states 
are different symmetries of the defects and the rates be- 


Experiment - Interstitials [xlO 2 ] 


-1.5014 

0.0000 

0.2827 

0.4473 

0.7714 

0.0040 

-2.1958 

2.1918 

0.0000 

0.0000 

0.0218 

0.5257 

-0.8295 

0.2820 

0.0000 

0.1171 

0.0000 

0.2219 

-0.4665 

0.1275 

0.3825 

0.0000 

0.0000 

0.4346 

-0.8172 

Simulation - Interstitials 

[xlO 2 ] 


-2.0873 

0.0000 

0.5688 

0.5834 

0.9351 

0.1637 

-1.9832 

1.6178 

0.2016 

0.0000 

0.0231 

0.1631 

-0.3263 

0.1401 

0.0000 

0.0578 

0.0000 

0.1093 

-0.3682 

0.2012 

0.2939 

0.0000 

0.0543 

0.4219 

-0.7701 

Experiment - Vacancies 

[xlO 2 ] 


-1.4511 

0.1469 

1.0814 

0.2227 

0.0000 

0.0694 

-0.9595 

0.5917 

0.2984 

0.0000 

0.1820 

0.1521 

-0.3341 

0.0000 

0.0000 

0.0509 

0.2643 

0.0000 

-0.5841 

0.2689 

0.0000 

0.0000 

0.0000 

1.0602 

-1.0603 

Simulation - Vacancies [ 

xlO 2 ] 


-1.7116 

0.2209 

0.9771 

0.2794 

0.2342 

0.0000 

-0.4375 

0.2625 

0.1734 

0.0015 

0.0874 

0.0812 

-0.1686 

0.0000 

0.0000 

0.0000 

0.1483 

0.0000 

-0.3230 

0.1747 

0.1621 

0.0000 

0.3248 

0.7269 

-1.2137 


TABLE II: Rate constant matrix Kij for F = 140 obtained 
from experiment and simulation using the optimization of Eq. 

El 


In summary we have shown that geometrical defects 
like interstitials and vacancies appear as different topo¬ 
logical configurations mainly constructed of two, three, 
or four dislocations where the burger-vector cancels. We 
characterize the configurations by the symmetry of 2D 
point groups and show that the symmetries for inter¬ 
stitials and vacancies are not equivalent. The relative 
equilibrium probabilities of defects vary as a function of 
the temperature. In the low temperature limit, the prob¬ 
abilities of 1 D I and 1 symmetries are largest for in¬ 
terstitials, while vacancies are predominantly in the v C\ 
symmetry. This completely different temperature depen¬ 
dence of vacancies and interstitials is not dominated by 
their energy, which is almost degenerated but by the en¬ 
tropy. The entropic and energetic contributions can be 


tween different states are determined from time corre¬ 
lation functions, which we measure in experiment and 
computer simulations. 

This work presents a detailed study on the defect ener¬ 
getics and dynamics of point defects in two dimensional 
materials. We hope that this motivates experiments in 
other two dimensional systems. A particular future ques¬ 
tion which becomes accessible, e.g. in ultracold dipolar 
quantum gases [30. and graphene [31] is the role quantum 
fluctuations in the defect dynamics. 
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